##########################################################################
##########################################################################
##### Issues vs. Affect: How Do Elite and Mass Polarization Compare? #####
##########################################################################
##########################################################################

library(lattice)
library(latticeExtra)
library(MASS)
library(foreign)
library(gridExtra)
library(car)
library(TeachingDemos)

#################################################################################################

#####
## Distributions of issue attitude- 
## and thermometer-based ideal points
## by party identification and political
## stratum in 2004
#####

ideodata <- read.dta("Combined, Short.dta")
head(ideodata)

ideodata2004 <- subset(ideodata, ideodata$year=="2004" & ideodata$rep!="NA")

ideodata2004$sample2[ideodata2004$sample==2] <- 1
ideodata2004$sample2[ideodata2004$sample==1] <- 2

plot1 <- densityplot(~issuescore | samplealt, 
                     data = ideodata2004,
                     #auto.key = TRUE,
                     aspect = .75,
                     ylab = "Density Estimate",
                     xlab = "Issue-Based Ideal Points",
                     main = "A. Ideological Polarization",
                     #ylim = c(0, 1.5),
                     groups = repalt,
                     scond = ideodata2004$sample2,
                     plot.points = FALSE,
                     layout = c(1, 2),
                     col = c("blue", "red"),
                     panel = panel.superpose,
                     panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
                       pnl = panel.number()
                       panel.densityplot(x, col=...)
                       panel.abline(v = -.948779[pnl==1], lty = 2, 
                                    lwd = .5, col = c("blue")[pnl==1])
                       panel.abline(v = .8013903[pnl==1], lty = 2, 
                                    lwd = .5, col = c("red")[pnl==1])
                       panel.abline(v = -.2340741[pnl==2], lty = 2, 
                                    lwd = .5, col = c("blue")[pnl==2])
                       panel.abline(v = .4156335[pnl==2], lty = 2, 
                                    lwd = .5, col = c("red")[pnl==2])
                     }
)

plot2 <- densityplot(~feelingscore | as.factor(samplealt), 
                     data = ideodata2004,
                     #auto.key = TRUE,
                     aspect = .75,
                     ylab = "Density Estimate",
                     xlab = "Thermometer-Based Ideal Points",
                     main = "B. Affective Polarization",
                     #ylim = c(0, 1.5),
                     groups = repalt,
                     scond = ideodata2004$sample2,
                     plot.points = FALSE,
                     layout = c(1, 2),
                     col = c("blue", "red"),
                     panel = panel.superpose,
                     panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
                       pnl = panel.number()
                       panel.densityplot(x, col=...)
                       panel.abline(v = -1.036329[pnl==1], lty = 2, 
                                    lwd = .5, col = c("blue")[pnl==1])
                       panel.abline(v = 1.117478[pnl==1], lty = 2, 
                                    lwd = .5, col = c("red")[pnl==1])
                       panel.abline(v = -.4685708[pnl==2], lty = 2, 
                                    lwd = .5, col = c("blue")[pnl==2])
                       panel.abline(v = .7483635[pnl==2], lty = 2, 
                                    lwd = .5, col = c("red")[pnl==2])
                     }
)          

print(plot1, position = c(0, 0, .5, 1), more = T)
print(plot2, position = c(.5, 0, 1, 1), more = F)

#################################################################################################

#####
## Average difference between party ideal
## point means over time, by political stratum
#####

meandata2 <- read.csv("Mean Differences, Short.csv", header = TRUE)
head(meandata2)

meandata2$cond[meandata2$item=="Issue-Based"] <- 1
meandata2$cond[meandata2$item=="Thermometer-Based"] <- 2

xyplot(estimate ~ year | item,
       data = na.omit(meandata2),
       aspect = .75,
       ylab = "Party Mean Difference",
       xlab = "Presidential Election Year",
       main = "A. Difference in Means",
       group = sample,
       scond = meandata2$cond,
       type = "b",
       layout = c(2, 1),
       ylim = c(-0.1,2.7),
       col = "black",
       pch = c(1, 16),
       #key=list(columns=1,points=list(pch=c(16,1)), 
       #        text=list(c("Masses", "Elites"))),
       par.settings = list(superpose.symbol = list(pch=c(16, 1), col="black")),           
       panel = panel.superpose,
       panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
         pnl = panel.number()
         panel.xyplot(x, y, pch=...)
         panel.segments(x0=meandata2$year[scond==pnl], y0=meandata2$lower[scond==pnl],
                        x1=meandata2$year[scond==pnl], y1=meandata2$upper[scond==pnl], col = "black")
       }
)

#################################################################################################

#####
## Difference in party medians,
## by political stratum
#####

mediandata <- read.csv("Median Differences, Short.csv", header = TRUE)
head(mediandata)

mediandata$cond[mediandata$item=="Issue-Based"] <- 1
mediandata$cond[mediandata$item=="Thermometer-Based"] <- 2

xyplot(estimate ~ year | item,
       data = na.omit(mediandata),
       aspect = .75,
       ylab = "Party Median Difference",
       xlab = "Presidential Election Year",
       main = "B. Difference in Medians",
       group = sample,
       scond = mediandata$cond,
       type = "b",
       layout = c(2, 1),
       ylim = c(-.1, 2.8),
       col = "black",
       pch = c(1, 16),
       #key=list(columns=1,points=list(pch=c(16,1)), 
       #        text=list(c("Masses", "Elites"))),
       par.settings = list(superpose.symbol = list(pch=c(16, 1), col="black")),           
       panel = panel.superpose,
       panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
         pnl = panel.number()
         panel.xyplot(x, y, pch=...)
         panel.segments(x0=mediandata$year[scond==pnl], y0=mediandata$lower[scond==pnl],
                        x1=mediandata$year[scond==pnl], y1=mediandata$upper[scond==pnl], col = "black")
       }
)

#################################################################################################

#####
## Average difference between party ideal
## point standard deviations over time,
## by political stratum
#####

attach(ideodata)

sqrt(sigma.test(ideodata[ which(sample==1 & year==1972),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1976),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1980),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1984),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1988),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1992),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1996),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2000),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2004),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2008),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2012),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2016),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])

sqrt(sigma.test(ideodata[ which(sample==2 & year==1972),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==1980),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==1984),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==1988),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==1992),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==2000),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==2004),]$feelingscore, 
                conf.level=0.95)$conf.int[1:2])

sqrt(sigma.test(ideodata[ which(sample==1 & year==1972),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1976),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1980),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1984),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1988),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1992),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==1996),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2000),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2004),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2008),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2012),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==1 & year==2016),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])

sqrt(sigma.test(ideodata[ which(sample==2 & year==1980),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==1984),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==1988),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==1992),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==2000),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(ideodata[ which(sample==2 & year==2004),]$issuescore, 
                conf.level=0.95)$conf.int[1:2])

stdevdata2 <- read.csv("Standard Deviations, Short.csv", header = TRUE)
head(stdevdata2)

stdevdata2$cond[stdevdata2$item=="Issue-Based"] <- 1
stdevdata2$cond[stdevdata2$item=="Thermometer-Based"] <- 2

xyplot(estimate ~ year | item,
       data = na.omit(stdevdata2),
       aspect = .75,
       ylab = "Standard Deviation",
       xlab = "Presidential Election Year",
       main = "C. Standard Deviations",
       group = sample,
       scond = stdevdata2$cond,
       ylim = c(0.4,1.5),
       type = "b",
       col = "black",
       pch = c(1, 16),
       layout = c(2,1),
       #key=list(columns=1,points=list(pch=c(16,1)), 
       #        text=list(c("Masses", "Elites"))),
       par.settings = list(superpose.symbol = list(pch=c(16, 1), col="black")),           
       panel = panel.superpose,
       panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
         pnl = panel.number()
         panel.xyplot(x, y, pch=...)
         panel.segments(x0=stdevdata2$year[scond==pnl], y0=stdevdata2$lower[scond==pnl],
                        x1=stdevdata2$year[scond==pnl], y1=stdevdata2$upper[scond==pnl], col = "black")
       }
)

#################################################################################################

#####
## Note: the following code replicates
## figures from the Supplemental Appendix
#####

#################################################################################################

#####
## Absolute difference in feeling 
## thermometer scores for "opposite"
## partisan/ideological stimuli over
## time for masses and elites
#####

thermdata3 <- read.csv("Thermometer Differences, Short.csv", header = TRUE)
head(thermdata3)

thermdata3$cond[thermdata3$item=="Conservatives"] <- 1
thermdata3$cond[thermdata3$item=="Democratic Candidate"] <- 2
thermdata3$cond[thermdata3$item=="Democratic Party"] <- 3
thermdata3$cond[thermdata3$item=="Liberals"] <- 4
thermdata3$cond[thermdata3$item=="Republican Candidate"] <- 5
thermdata3$cond[thermdata3$item=="Republican Party"] <- 6

xyplot(estimate ~ year | as.factor(item), 
       data = na.omit(thermdata3),
       aspect = 0.75,
       ylab = "Mean Party Difference",
       xlab = "Presidential Election Year",
       ylim = c(0.2,2.2),
       main = "",
       col = "black",
       type = "b",
       group = sample,
       scond = thermdata3$cond,
       layout = c(3, 2),
       pch = c(1, 16),
       key=list(columns=1,points=list(pch=c(16, 1)), 
                text=list(c("Masses", "Elites"))),
       panel = panel.superpose,
       panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
         pnl = panel.number()
         panel.xyplot(x, y, pch=...)
         panel.segments(x0=thermdata3$year[scond==pnl], y0=thermdata3$lower[scond==pnl],
                        x1=thermdata3$year[scond==pnl], y1=thermdata3$upper[scond==pnl], col = "black")
       }
)

#################################################################################################

#####
## Average difference between Republican
## and Democratic self-identifiers attitudes
## about several issues over time for masses
## and elites
#####

issuedata2 <- read.csv("Issue Differences, Short.csv", header = TRUE)
head(issuedata2)

issuedata2$cond[issuedata2$item=="Abortion"] <- 1
issuedata2$cond[issuedata2$item=="Aid to Minorities"] <- 2
issuedata2$cond[issuedata2$item=="Defense Spending"] <- 3
issuedata2$cond[issuedata2$item=="Government Services"] <- 4
issuedata2$cond[issuedata2$item=="Health Insurance"] <- 5

issuedata2$yearalt <- factor(issuedata2$year)

xyplot(estimate ~ year | as.factor(item), 
       aspect = .75,
       data = na.omit(issuedata2),
       ylab = "Mean Party Difference",
       xlab = "Presidential Election Year",
       ylim = c(-0.2,2),
       main = "",
       col = "black",
       type = "b",
       group = sample,
       scond = issuedata2$cond,
       layout = c(3, 2),
       pch = c(1,16),
       #scales=list(x = list(labels = c("'88", "'92", "'00", "'04", "'12"))),
       key=list(columns=1,points=list(pch=c(16,1)), 
                text=list(c("Masses", "Elites")), corner=c(.9,.75)),
       panel = panel.superpose,
       panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
         pnl = panel.number()
         panel.xyplot(x, y, pch=...)
         panel.segments(x0=issuedata2$year[scond==pnl], y0=issuedata2$lower[scond==pnl],
                        x1=issuedata2$year[scond==pnl], y1=issuedata2$upper[scond==pnl], col = "black")
       }
)

#################################################################################################

#####
## Average difference between parties 
## using standardized index measure
#####

meandata2 <- read.csv("Mean Differences, Index.csv", header = TRUE)
head(meandata2)

meandata2$cond[meandata2$item=="Issue-Based"] <- 1
meandata2$cond[meandata2$item=="Thermometer-Based"] <- 2

xyplot(estimate ~ year | item,
       data = na.omit(meandata2),
       aspect = .75,
       ylab = "Mean Party Difference",
       xlab = "Presidential Election Year",
       group = sample,
       scond = meandata2$cond,
       type = "b",
       ylim = c(-.1, 2),
       layout = c(2, 1),
       col = "black",
       pch = c(1, 16),
       key=list(columns=1,points=list(pch=c(16,1)), 
                text=list(c("Masses", "Elites"))),
       par.settings = list(superpose.symbol = list(pch=c(16, 1), col="black")),           
       panel = panel.superpose,
       panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
         pnl = panel.number()
         panel.xyplot(x, y, pch=...)
         panel.segments(x0=meandata2$year[scond==pnl], y0=meandata2$lower[scond==pnl],
                        x1=meandata2$year[scond==pnl], y1=meandata2$upper[scond==pnl], col = "black")
       }
)

#################################################################################################

#####
## Index standard deviations over time, 
## by political stratum
#####

attach(ideodata)

sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1972),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1976),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1980),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1984),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1988),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1992),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1996),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2000),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2004),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2008),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2012),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2016),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])

sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1972),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1980),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1984),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1988),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1992),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==2000),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==2004),]$feelingindex), 
                conf.level=0.95)$conf.int[1:2])

sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1972),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1976),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1980),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1984),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1988),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1992),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==1996),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2000),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2004),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2008),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2012),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==1 & year==2016),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])

sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1980),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1984),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1988),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==1992),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==2000),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])
sqrt(sigma.test(na.omit(ideodata[ which(sample==2 & year==2004),]$issueindex), 
                conf.level=0.95)$conf.int[1:2])

stdevdata2 <- read.csv("Standard Deviations, Index.csv", header = TRUE)
head(stdevdata2)

stdevdata2$cond[stdevdata2$item=="Issue-Based"] <- 1
stdevdata2$cond[stdevdata2$item=="Thermometer-Based"] <- 2

xyplot(estimate ~ year | item,
       data = na.omit(stdevdata2),
       aspect = .75,
       ylab = "Standard Deviation",
       xlab = "Presidential Election Year",
       group = sample,
       scond = stdevdata2$cond,
       type = "b",
       col = "black",
       ylim = c(0.49, 1.15),
       pch = c(1, 16),
       layout = c(2,1),
       key=list(columns=1,points=list(pch=c(16,1)), 
                text=list(c("Masses", "Elites"))),
       par.settings = list(superpose.symbol = list(pch=c(16, 1), col="black")),           
       panel = panel.superpose,
       panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
         pnl = panel.number()
         panel.xyplot(x, y, pch=...)
         panel.segments(x0=stdevdata2$year[scond==pnl], y0=stdevdata2$lower[scond==pnl],
                        x1=stdevdata2$year[scond==pnl], y1=stdevdata2$upper[scond==pnl], col = "black")
       }
)

#################################################################################################

#####
## Party medians for index measure,
## by political stratum
#####

mediandata <- read.csv("Median Differences, Index.csv", header = TRUE)
head(mediandata)

mediandata$cond[mediandata$item=="Issue-Based"] <- 1
mediandata$cond[mediandata$item=="Thermometer-Based"] <- 2

xyplot(estimate ~ year | item,
       data = na.omit(mediandata),
       aspect = .75,
       ylab = "Party Median Difference",
       xlab = "Presidential Election Year",
       group = sample,
       scond = mediandata$cond,
       type = "b",
       layout = c(2, 1),
       ylim = c(-.1, 2.2),
       col = "black",
       pch = c(1, 16),
       key=list(columns=1,points=list(pch=c(16,1)), 
                text=list(c("Masses", "Elites"))),
       par.settings = list(superpose.symbol = list(pch=c(16, 1), col="black")),           
       panel = panel.superpose,
       panel.groups = function(x, y, subscripts, groups, col, lty, scond, ...) {
         pnl = panel.number()
         panel.xyplot(x, y, pch=...)
         panel.segments(x0=mediandata$year[scond==pnl], y0=mediandata$lower[scond==pnl],
                        x1=mediandata$year[scond==pnl], y1=mediandata$upper[scond==pnl], col = "black")
       }
)

#################################################################################################

######
## Plot joint scaling model discrimination
## and difficulty parameters
#####

# Discrimination, feeling thermometers
disctherms <- read.csv("Therms Discrimination, Short.csv", header = TRUE)
head(disctherms)

dotplot(as.factor(item) ~ joint, 
        data = disctherms,
        aspect = 1.5,
        xlim = c(-3.5, 4.5),
        xlab = "Discrimination Parameter",
        ylab = "",
        #main = "A. Feeling Thermometers",
        panel=function(...){
          panel.dotplot(..., lty = 2)
          panel.xyplot(disctherms$joint, disctherms$item, col = "black", pch = 16)
          panel.segments(disctherms$lowerj, as.numeric(disctherms$item),
                         disctherms$upperj, as.numeric(disctherms$item), lty = 1, col = "black")
          panel.abline(v=0, lty=2)
        }
)

# Discrimination, issue attitudes
discissues <- read.csv("Issues Discrimination, Short.csv", header = TRUE)
head(discissues)

dotplot(as.factor(item) ~ joint, data = discissues,
        aspect = 1.5,
        xlim = c(-3.5, 4.5),
        xlab = "Discrimination Parameter",
        ylab = "",
        #main = "B. Issue Attitudes",
        col = 1,
        lyt = 2,
        panel=function(...){
          panel.dotplot(..., lty = 2)
          panel.xyplot(discissues$joint, discissues$item, col = "black", pch = 16)
          panel.segments(discissues$lowerj, as.numeric(discissues$item),
                         discissues$upperj, as.numeric(discissues$item), lty = 1, col = "black")
          panel.abline(v=0, lty=2)
        }
)


# Difficulty, feeling thermometers
difftherms <- read.csv("Therms Difficulty, Short.csv", header = TRUE)
head(difftherms)

dotplot(cutpoint ~ diff, 
        data = difftherms,
        aspect = 1.5,
        #xlim = c(.4, .6),
        xlab = "Difficulty Parameter",
        ylab = "",
        col = 1,
        lyt = 2,
        scales=list(
          y=list(labels = rownames(difftherms$cutpoint))),
        panel=function(x, y){
          panel.dotplot(x, y, col = "black", lty = 2)
          #panel.segments(diff.data$lower, as.numeric(y),
          #              + diff.data$upper, as.numeric(y), lty = 1, col = "black")
          #panel.abline(v=0, lty=2)
        }
)


# Difficulty, issue attitudes
diffissues <- read.csv("Issues Difficulty, Short.csv", header = TRUE)
head(diffissues)

dotplot(cutpoint ~ diff, data = diffissues,
        aspect = 1.5,
        #xlim = c(.4, .6),
        xlab = "Difficulty Parameter",
        ylab = "",
        col = 1,
        lyt = 2,
        scales=list(
          y=list(labels = rownames(diffissues$cutpoint))),
        panel=function(x, y){
          panel.dotplot(x, y, col = "black", lty = 2)
          #panel.segments(diff.data$lower, as.numeric(y),
          #              + diff.data$upper, as.numeric(y), lty = 1, col = "black")
          #panel.abline(v=0, lty=2)
        }
)

#################################################################################################

######
## Provide robustness checks for
## joint scaling models per advice
## of Jessee (2016)
#####

ideodata <- read.dta("Combined, Short.dta")
head(ideodata)

#####
## First, compare ideal points
## across estimations, by item type
#####

massonly <- subset(ideodata, ideodata$sample == 1, )
eliteonly <- subset(ideodata, ideodata$sample == 2, )

cor(massonly$feelingsANES, massonly$feelingscore)
cor(massonly$issuesANES, massonly$issuescore)

xyplot(feelingsANES ~ feelingscore, data = massonly,
       aspect = 1,
       ylab = "Mass Estimation",
       xlab = "Joint Estimation",
       main = "A. Feeling Thermometers",
       col = "light grey",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")
       }
)

xyplot(issuesANES ~ issuescore, data = massonly,
       aspect = 1,
       ylab = "Mass Estimation",
       xlab = "Joint Estimation",
       main = "B. Issue Attitudes",
       col = "light grey",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")
       }
)


cor(eliteonly$feelingsCDS, eliteonly$feelingscore)
cor(eliteonly$issuesCDS, eliteonly$issuescore)

xyplot(feelingsCDS ~ feelingscore, data = eliteonly,
       aspect = 1,
       ylab = "Elite Estimation",
       xlab = "Joint Estimation",
       main = "C. Feeling Thermometers",
       col = "light grey",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")         
       }
)

xyplot(issuesCDS ~ issuescore, data = eliteonly,
       aspect = 1,
       ylab = "Elite Estimation",
       xlab = "Joint Estimation",
       main = "D. Issue Attitudes",
       col = "light grey",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")         
       }
)

#####
## Next, compare discrimination 
## parameters across estimations,
## by item type
#####

# Feeling thermometers
disctherms <- read.csv("Therms Discrimination, Short.csv", header = TRUE)
head(disctherms)

cor(disctherms$joint, disctherms$mass)

xyplot(joint ~ -1*mass, data = disctherms,
       aspect = 1,
       ylab = "Joint Estimation",
       xlab = "Mass Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         #panel.loess(x, y, lty = 1, col = "red")
       }
)


cor(disctherms$joint, disctherms$elite)

xyplot(joint ~ -1*elite, data = disctherms,
       aspect = 1,
       ylab = "Joint Estimation",
       xlab = "Elite Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         #panel.loess(x, y, lty = 1, col = "red")         
       }
)


cor(disctherms$mass, disctherms$elite)

xyplot(mass ~ elite, data = disctherms,
       aspect = 1,
       ylab = "Mass Estimation",
       xlab = "Elite Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         #panel.loess(x, y, lty = 1, col = "red")         
       }
)

# Issue attitudes
discissues <- read.csv("Issues Discrimination, Short.csv", header = TRUE)
head(discissues)

cor(discissues$joint, discissues$mass)

xyplot(joint ~ mass, data = discissues,
       aspect = 1,
       ylab = "Joint Estimation",
       xlab = "Mass Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         #panel.loess(x, y, lty = 1, col = "red")
       }
)


cor(discissues$joint, discissues$elite)

xyplot(joint ~ elite, data = discissues,
       aspect = 1,
       ylab = "Joint Estimation",
       xlab = "Elite Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         #panel.loess(x, y, lty = 1, col = "red")         
       }
)


cor(discissues$mass, discissues$elite)

xyplot(mass ~ elite, data = discissues,
       aspect = 1,
       ylab = "Mass Estimation",
       xlab = "Elite Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         #panel.loess(x, y, lty = 1, col = "red")         
       }
)

#####
## Finally, compare difficulty 
## parameters across estimations,
## by item type
#####

# Feeling thermometers
difftherms <- read.csv("Therms Difficulty, Short.csv", header = TRUE)
head(difftherms)

cor(difftherms$diff, difftherms$mass)

xyplot(diff ~ -1*mass, data = difftherms,
       aspect = 1,
       ylab = "Joint Estimation",
       xlab = "Mass Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")
       }
)


cor(difftherms$diff, difftherms$elite)

xyplot(diff ~ -1*elite, data = difftherms,
       aspect = 1,
       ylab = "Joint Estimation",
       xlab = "Elite Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")         
       }
)


cor(difftherms$mass, difftherms$elite)

xyplot(mass ~ elite, data = difftherms,
       aspect = 1,
       ylab = "Mass Estimation",
       xlab = "Elite Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")         
       }
)

# Issue attitudes
diffissues <- read.csv("Issues Difficulty, Short.csv", header = TRUE)
head(diffissues)

cor(diffissues$diff, diffissues$mass)

xyplot(diff ~ mass, data = diffissues,
       aspect = 1,
       ylab = "Joint Estimation",
       xlab = "Mass Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")
       }
)


cor(diffissues$diff, diffissues$elite)

xyplot(diff ~ elite, data = diffissues,
       aspect = 1,
       ylab = "Joint Estimation",
       xlab = "Elite Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")         
       }
)


cor(diffissues$mass, diffissues$elite)

xyplot(mass ~ elite, data = diffissues,
       aspect = 1,
       ylab = "Mass Estimation",
       xlab = "Elite Estimation",
       main = "",
       col = "black",
       panel=function(x, y, ...){
         panel.xyplot(x, y, ...)
         panel.lmline(x, y, lty = 1, col = "black")
         panel.loess(x, y, lty = 1, col = "red")         
       }
)